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ABSTRACT 

Atmospheres having a significant radiative support are shown to be intrinsically unstable at luminosi- 
ties above a critical fraction Pdit ~ 0.5 — 0.85 of the Eddington limit, with the exact value depending on 
the boundary conditions. Two different types of absolute radiation-hydrodynamic instabilities of acous- 
tic waves are found to take place even in the electron scattering dominated limit. Both instabilities grow 
over dynamical time scales and both operate on non radial modes. One is stationary and arises only 
after the effects of the boundary conditions are taken into account, while the second is a propagating 
wave and is insensitive to the boundary conditions. Although a significant wind can be generated by 
these instabilities even below the classical Eddington luminosity limit, quasi-stable configurations can 
exist beyond the Eddington limit due to the generally reduced effective opacity. 

The study is done using a rigorous numerical linear analysis of a gray plane parallel atmosphere under 
the Eddington approximation. We also present more simplified analytical explanations. 

Subject headings: Radiative transfer — hydrodynamics — instabilities — stars: atmospheres — stars: 
variables: other — stars: oscillations 



1. INTRODUCTION 

Luminous objects shining close to the Eddington limit 
are of special interest since they occur in a wide range of 
astrophysical contexts. As such, they have been a popu- 
lar topic for observational and theoretical investigations. 
Luminosities in the range of the Eddington limit appear 
in the very luminous stars populating the top of the HR 
diagram, in novae during the critical phases of the erup- 
tion and also in accreting objects of both galactic and 
extra-galactic scales. These very luminous objects usu- 
ally display a great diversity of phenomena ranging from 
erratic variability to spatial structures and winds. Con- 
sequently, the theoretical modeling requires the incorpo- 
ration of hydrodynamics, radiation and sometimes even 
magnetic fields. Undoubtedly, much of the interesting be- 
havior stems from the large, dynamically important lumi- 
nosity of these objects - luminosities that are close to the 
Eddington limit. 

The classical Eddington limit is the upper limit possible 
to the luminosity of an object in hydrostatic balance (Ed- 
dington 1926). It is the luminosity for which the radiative 
force upwards balances the gravitational force downwards 
assuming the opacity of the (fully ionized) gas attains its 
minimal possible value, that is, the Thomson scattering 
opacity for non relativistic electrons. For a homogeneous 
static configuration, it is the upper luminosity limit be- 
cause above it there is no consistent solution to the hydro- 
static equations. An atmosphere that violates this upper 
limit will then simply blow itself apart - an optically very 
thick wind is then unavoidable and the object will evap- 
orate on a relatively short time scale. When the opac- 
ity is larger than the Thomson scattering one, as is often 
the case, one can define a modified Eddington limit that 
is smaller than the 'classical' value (e.g., Humphreys & 



Davidson 1984, Lamers 1986, and Appenzeller 1986). The 
modified Eddington luminosity was used to explain the 
upper luminosity limit observed for stars and the variabil- 
ity that these stars exhibit. The meaning of the modified 
Eddington luminosity should however be taken cautiously, 
because as one approaches the modified Eddington limit, 
the increased scale height and reduced density imply that 
in any fully ionized atmosphere, the opacity law is reduced 
to that of Thomson scattering (Humphreys & Davidson 
1994). Namely, as the modified limit is approached, it it- 
self approaches the classical limit such that the modified 
Eddington limit is probably not an absolute luminosity 
limit. 

On top of the interesting questions on the nature of at- 
mospheres when their luminosity approaches the Edding- 
ton luminosity, exists the fundamental question: Ts the 
Eddington luminosity the maximal luminosity under all 
conditions?'. Part of the answer was already answered in 
the negative by Shaviv (1998). Here we elaborate and ex- 
tend the analysis and demonstrate that the classical Ed- 
dington limit is the limit under a very restrictive set of 
assumptions. Still we would like to emphasize that a sig- 
nificant factor in turning the Eddington limit into a major 
tool in setting astrophysical limits is the simplicity of the 
classical expression, namely Ce = 47rcGA//CTThomson- 

There are several reasons why the study of radiative in- 
stabilities is important. Most, if not all of the objects shin- 
ing close to the Eddington limit show variability. Novae 
for example, not only show very high (occasionally super- 
Eddington) luminosities over extended periods of time, 
they also show clumpy ejected material (e.g., the high res- 
olution HST images, Shara et al. 1997) which could be 
explained using radiatively driven instabilities. Another 
example are the super-Giants populating the top of the 
HR diagram. Quite a few stars with M ~ lOOM© and 
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L ~ lO^Z/0 exist across a wide range of spectral types and 
exhibit various degrees of variability. For example, several 
of the numerous 03 type stars at 30 Doradus are so bright 
that their luminosity is clearly very close to their Edding- 
ton limit (e.g., Massey & Hunter 1998). Less blue (and less 
young) are the Luminous Blue Variables (LBVs) termed 
so because of their often eruptive behavior in which their 
luminosity and rate of mass ejection can increase while 
their effective temperature decreases (see for example the 
review by Humphreys & Davidson 1994). When these ob- 
jects are in their eruptive state, at least some of them 
clearly surpass the Eddington limit for significant lengths 
of time. For example, 77 Car in its eruption between 1830 
and 1860 released roughly 10^^ ergs of light over a period 
much longer than its dynamical time scale (e.g., van Gon- 
deren & The 1985), yielding an average luminosity that 
is more than an order of magnitude higher than the Ed- 
dington limit (provided that the mass of 77 Car is of the 
order of 1{){)Mq). The variability of these objects which is 
probably a consequence of the high luminosity, is therefore 
interesting to understand. 

An additional type of interesting objects are 'very' and 
'super massive objects' (VMO's and SMO's) that might 
have existed as massive pop HI stars in the young universe. 
These massive stars are supposed to shine at luminosities 
very near the Eddington limit with the more massive of 
them at virtually the Eddington limit itself. Since the 
mctalicity of these objcc;ts is extremely small, the analysis 
carried out here, which assumes that the opacity is dom- 
inated by Thomson scattering, is particularly relevant to 
these objects. 

Super-critical accretion disks are an altogether different 
type of relevant luminous objects. Irrespective on whether 
they advect most of the excess accreted mass and energy 
above the Eddington rate into the black hole or whether 
the excess mass is blown away, each element of these disks 
should shine very close to local effective Eddington limit 
as these disks are always radiation pressure dominated. 

Various works have shown that atmospheres shining 
close to the Eddington limit exhibit a wide range of in- 
stabilities as well as other phenomena such as the gen- 
eration of strong winds. For example, a large radiation 
pressure contribution reduces the adiabatic index. The 
relatively 'loose' structure becomes unstable against con- 
vection when the adiabatic index decreases due to ioniza- 
tion (e.g., Ledoux 1965, Stothers & Chin 1993, and Wa- 
genhuber & Weiss 1994). 

Joss, Salpeter & Ostriker (1973) analyzed the appear- 
ance of convection at high stellar luminosities assuming 
that the condition for convective instability is the classical 
Schwarzchild condition, namely the temperature gradient 
must be higher than the adiabatic one. The main concern 
of Joss et al. was the energy flux and its relation to the 
Eddington luminosity. Here we are predominately inter- 
ested in the conditions for the stability under high radia- 
tive fluxes. The main result of Joss et al. was that as the 
energy flux through an atmosphere is increased, convection 
arises before the Eddington limit is reached. Moreover, the 
convective flux will always grow enough such that the re- 
maining radiative flux is less than Eddington, 05 long as 
convection is efficient. Namely, the interior of stars which 
are dense, will never witness super Eddington fltixes. It is 



only near the surface where super Eddington fluxes can be 
reached. Since the optical depth below which convection 
becomes inefficient is relatively high, the wind of a super 
Eddington object should be very thick and the mass loss 
should be enormous if nothing happens above the convec- 
tive zone (cf Shaviv 2000a). 

When the radiative force term becomes as important as 
the gravity or the gas pressure gradient terms, it can be 
a source of mechanical types of instabilities. For exam- 
ple, Glatzcl & Kiriakidis (1993) have shown that strange 
modes (which are 'opacity' modified acoustic waves) suffer 
an instability in which mechanical work by the radiation is 
pumped into the waves. This instability arises from the in- 
teraction of the acoustic waves with the non local nature of 
the diffusive radiation equations. Another instability that 
requires special opacity laws is the K-mechanism. Unlike 
the s-mode instability, however, it grows on time scales 
much slower than the dynamical time scale, it docs not 
involve work done by the radiation, and exists only when 
the system is neither fully adiabatic nor fully in the oppo- 
site (NAR) limit. However, it can exist when the system 
is far from Eddington, because it involves the transfer of 
heat and not work by the radiation. 

A related instability arises when the radiation field 
transfers energy into acoustic waves due to high radiation 
pressures, as was shown by Hearn (1972, 1973) and elabo- 
rated by Berthomieu et al. (1976), Spiegel (1976), Carlberg 
(1980) and Asplund (1998). The instability arises when 
the opacity increases under compression and the change in 
the opacity is synchronized with the radiative force. This 
instability is not a real instability (aka 'an absolute insta- 
bility') in the sense that an infinitesimal perturbation will 
not grow to become nonlinear, however, as it is a secular 
drift instability (aka 'a convective instability'), it can take 
finite amplitude waves (excited by convection for example) 
and increase their amplitude by many e-folds during one 
crossing time of the system. (A counter propagating wave 
will be dissipated as fast, thus adding up to no net am- 
plification of a standing wave unless the conditions for the 
strange mode instability are satisfied, e.g., Glatzel 1994, 
Papaloizou et al. 1997). An additional term arises when 
the Doppler variations of the radiative absorption (arising 
from wave motion) are taken into account (Shaviv 1999b). 
A third term is the driving through line absorption. This 
term which utilizes the Doppler effect as well, is often very 
important in wind acceleration from hot stars (Owocki & 
Rybicki 1984, 1985, 1986, 1991, Rybicki, Owocki & Cas- 
tor 1990, Feldmeier et al. 1997). The modes in systems 
with this type of simple drift instabilities cannot be used 
however to compose unstable standing waves in a simple 
Thomson scattering dominated atmosphere. 

The problem of the behavior of luminous atmospheres 
has also been compared to fluidized beds (Prendergast & 
Spiegel 1973). Intuitively, the kind of bubbling observed 
in fluidized beds should also manifest itself in radiatively 
supported atmospheres. Additional effects arise with the 
introduction of strong magnetic fields. Specifically, when 
strong magnetic fields are present such that hydrodynamic 
motion perpendicular to the magnetic field lines is frozen, 
a linear analysis shows that 'photon bubbles' arise at a 
high luminosity, facilitating the transfer of energy (Arons 
1992). The original context of accretion onto magnetized 
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neutron stars was also used to describe accretion disks as 
well (Gammie 1998). Magnetic fields are neglected in the 
present work. As we shall see, a plethora of instabilities 
can already manifest themselves when the magnetic fields 
are dynamically not important. 

In a recent work by Spiegel & Tao (1999), it was shown 
using a global non-radial linear mode analysis, that at- 
mospheres can have unstable standing modes (as opposed 
to the secular drift instability, or convective instability) 
even in the case of pure Thomson scattciring (as opposed 
to strange mode instability). These unstable modes were 
shown to be horizontally propagating waves that amplify 
over a time scale much longer than the dynamical times 
scale of the system. We have tried to recover this finding 
but could not succeed. We report in the Appendix how 
such an instability could artificially arise. Moreover, the 
two other soon be described instabilities are dynamically 
much more important. 

In addition to the intrinsic variability and the genera- 
tion of winds, perhaps the most interesting consequence 
of the instabilities is the change induced in the Edding- 
ton limit. Shaviv (1998) has shown that in the presence 
of inhomogcneities, the spatially averaged flux through an 
atmosphere can increase without increasing the time av- 
eraged force on it. This arises because of the tendency 
of radiation to find the paths of least resistance when es- 
caping through an atmosphere. Under some opacity laws, 
this will translate into having more radiation in regions 
where there is less mass and a smaller average force will 
be exerted. Thus, the luminosity can exceed the Edding- 
ton luminosity while the average force on a parcel of mass 
can remain less than the critical one. 

In summary, the classical Eddington limit is restricted to 
complete homogeneous systems and is frequently not valid 
in inhomogeneous ones. A fundamental question of how 
such inhomogcneities are formed and what form they ac- 
cept is discussed only briefly in Shaviv (1998). Here we ex- 
pand on this question with emphasis on Thomson scatter- 
ing dominated atmospheres ('Thomson atmospheres' here- 
after). 

Evidently, a picture which relates all the different insta- 
bilities together and gives sound physical reasoning to their 
origin is mostly absent even though it is of paramount im- 
portance to the understanding of very luminous objects. It 
is therefore, our goal in this paper to address the following 
points: 

1. We wish to unravel all the absolute instabilities that 
Thomson atmospheres (which are the simplest of all) 
exhibit as they approach the Eddington limit. This is 
done using a rigorous linear analysis (that is similar to 
but not the same as that of Spiegel & Tao 1999) which 
assumes nothing save a gray slab atmosphere under the 
Eddington approximation. 

2. We wish to understand the physical basis behind the 
instabilities found and how systems with these insta- 
bilities behave. 

Two instabilities will be found and described. Addi- 
tional instabilities arise when the effects of a varying opac- 
ity are taken into account. As the Eddington limit is 
approached, these instabilities generally become less im- 
portant as the opacity becomes mostly that of scattering. 



However, a future 'tabulation' of all the instabilities of lu- 
minous atmospheres should include them as well. 

A few clear extensions, such as a complete nonlinear 
analysis using hydrodynamic simulations or a systematic 
study of the parameter space, will be analyzed in forthcom- 
ing publications. These will also include a more elaborate 
global analysis of the non-trivial opacity law triggered in- 
stabilities. 

The paper is organized in the following way: The gov- 
erning radiative hydrodynamic equations are presented in 
§2. To clarify the physical picture, we discuss in §3 the pa- 
rameter space of the problem and show where and when 
the various limits appear. We then proceed in §4 to ana- 
lyze the problem in a full global analysis assuming only a 
gray plane parallel atmosphere, calculate the cigcnmodes, 
and find the characteristics of the instabilities found. We 
continue in §5 to describe the first instability and in §6 the 
second instability. We end with a discussion and summary 
of conclusions. An appendix is given to relate the results 
found here to those found by Spiegel & Tao (1999). 

2. THE RADIATIVE HYDRODYNAMIC EQUATIONS AND 
NOTATION 

We begin by laying out the equations used in this work. 
The equations describing the gas are the continuity, mo- 
mentum and energy equations: 

|^ + V-(pu) = (1) 

p— = ~Vp-pgz + p^F 2 
JJt c 

with standard notation: p, p and u are the gas (mat- 
ter) pressure, density and velocity, x is the total opac- 
ity per unit mass, which is the sum of the scattering a 
and the absorption k opacities. Throughout this paper, 
opacity coefficients without subscripts denote opacity per 
unit mass while a subscript v denotes opacity per unit vol- 
ume (with dimensions of length~^). S is the equilibrium 
radiation energy density calculated from the gas temper- 
ature S = AasBT^/c where usb is the Stefan Boltzmann 
constant. E is the radiation energy density and F is the 
radiation energy flux. The radiation energy and flux are 
those measured in the material rest frame, namely, the 
'Lagrangian' quantities. D/Dt = d/dt + u • V is the co- 
moving derivative. We distinguish between the adiabatic 
speed of sound = 7p/p and the isothermal speed of 
sound = p/p. Gravity acts downward along the z di- 
rection. 

In view of the above assumptions, we describe the ra- 
diation field using the Eddington approximation. Its 
'Lagrangian' equations, valid to 0{v/c), are (Mihalas & 
Weibel Mihalas 1984): 

DF 

— + -\/E = -p{K + a)cF. (5) 

The second term in the first equation is the 'compression' 
term which corrects the stationary frame radiation quan- 
tities into the material rest frame. This 'Doppler' term 
describes radiation drag. 
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The Eddington approximation is a reasonable assump- 
tion since we are interested in regions where oo > t > 1. 
Even at the photosphere where the approximation is least 
accurate, it yields reasonable quantitative results and the 
correct qualitative results (e.g., Mihalas & Weibel Miha- 
las 1984, pp. 357, 518). 

Once we apply a perturbation, quantities with an index 
'0' are the unperturbed values while an index '1' indicates 
the first order perturbation. The first order perturbation 
to the opacity, Xv,i, can be written as: 



Xv,l — Xv,0 



'- pi ^- Ti 

Xp—+XT' 



(6) 



Po "To 

with the following definitions for the logarithmic deriva- 
tives: 



Xp 



5 In Xv 



dlnp 



and xt 



9 In Xv 



d\nT 



(7) 



The pure Thomson scattering case corresponds to Xp = ^ 
and Xt = 0. However, the approximation we will use in 
this paper is 

^ k/o- < 1, (8) 

namely, the absorption is very small compared to the scat- 
tering, so that the derivatives arc to a good approximation 
those of pure Thomson. However, the absorption is not 
zero and the gas absorbs a small amount of radiation. This 
allows at times for the gas to equilibrate its temperature 
with the radiation, something that for a purely scattering 
atmosphere cannot take place. 

The last equation is the equation of state of a perfect 
gas relating the gas variables: 

Amu 

where A is the mean atomic weight of the gas and ks the 
Boltzmann constant. 

As the radiation field may not be in equilibrium with the 
matter, the expression for the opacity can be very compli- 
cated and in principle the opacity should be calculated 
consistently with the radiation field. In the present work 
we restrict the discussion to a gray opacity. 

The Fourier transform of a quantity x is denoted by x. 

3. DIMENSIONAL ANALYSIS AND CRITICAL PARAMETERS 

The driving forces on the system arc the terms that ap- 
pear on the r.h.s. of eq. (2). In equilibrium, the radiative 
and pressure forces balance the pull of gravity and con- 
sequently the r.h.s. vanishes. The relative importance of 
each term can be related to the ratio between the radiation 
pressure and the gas pressure: 

■Bo/3 

P 

where p = pgas • 

A second related parameter which determines the ratio 
between the radiative and gravitational forces is: 

Fx 



P _ Prad 



(10) 



r = 



eg 



(11) 



where F is the radiative fiux and x the opacity per unit 
mass, r is often called the Eddington parameter. If x 
is taken to be the Thomson scattering cross-section, then 



r = 1 corresponds to the Eddington limit. For other opac- 
ities (which are generally larger if the gas is ionized) F = 1 
corresponds to the Modified Eddington limit. 

The parameters and F are simply related to each other 
in Thomson atmospheres, for which the opacity per unit 
mass is constant. From the Eddington approximation we 
have that 



Fn 

(2 + 3r)^, 
c 



(12) 



where r is the optical depth for extinction. The decrease in 

the effective gravitational pull (1 — r)g leads to an increase 
in the density scale height of the atmosphere, namely. 



5(1 -r) 



(13) 



If the radiative force is constant and isothermality is as- 
sumed, then the puffed up scale height is constant with 
height. In the case of an atmosphere with an adiabatic 
gradient, ct should be replaced with c^, in which case the 
scale height is not constant with height any more. To 
demonstrate the relation between f3 and F we assume for 
simplicity an isothermal photosphere. In this approxima- 
tion, 



Poiz) = pq{z = 0) exp(-z/?p); 
Pq{z) = CtPo{z = 0) exp{-z/lp). 



(14) 



The optical depth from 2; = oo to a given point at a depth 
of z is 

POO 

T= PXdz = xlpPoiz = 0) exp{-z/lp) = xlpPo{z)/c%. 

J Z 

(15) 

For further reference we also have t = Xv^p- We thus find 
from eqs. (10) through (15) that 



^^Prad^(2 + 3T)^ 

■P 4:CTC^ 

r 

for T > 1. 



(2 + 3t) F 

3t 1-: 



(l-F) 



(16) 



Conversely, for an optically thick atmosphere, we have that 



(17) 



Atmospheres that have a small radiation to gas pressure 
ratio will therefore have a luminosity that is considerably 
smaller than the Eddington luminosity (given by F = 1). 
When the gas and radiation pressures equate, we have that 
F « 1/2. When the radiation pressure is much larger than 
the gas pressure, F approaches unity (such that 1 — F ^ 1). 
In general, when the opacity is a function of height or 
when the non-isothcrmalncss of the atmosphere is consid- 
ered, we will get a more complicated relation between F 
and (3, however, we will still find that /3 <C 1 corresponds 
to F <C 1, that /? « 1 corresponds to F ^ 1/2 and that 
/3 3> 1 corresponds to (1 — F) ^ 1. 

Additionally important parameters arise when compar- 
ing the different time scales found in the problem. The first 
time scale is the dynamic time scale defined here as the 
isothermal sound crossing time of a density scale height: 



= If. - cr 



(18) 
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From analysis of equations (4) and (5), a second time scale 
appears. To see this, we combine the two linearized ver- 
sions of the equations assuming the heat absorption is neg- 
ligible (i.e., in the adiabatic or isothermal limits) and ne- 
glecting the 'Doppler' term, to get: 

1 d^E XvdE I 

where Xv = per- In the more general case the coefficient of 
dE/dt is p{(j + k)/c. If the typical length scale is given by 
fc"^, a diffusion time scale can be defined as 



(20) 



when comparing the second and third terms. Specifically, 
for a wavelength of the order of the scale height of the 
atmospheres, we can define: 

3 Xvll 



Tdiff = Tdifi(A; = 2-1^1 p ) = 



{2-Kf 



(21) 



A dimensionless parameter is its ratio to the dynamical 
time scale, given by 

3 Xvll ct 



rdiff = 



Tdiff 
Tdyn 

3 



(27r)2 



(22) 



(2^ (t) ^""^^"^ = (2^ (?) 



where tq = Xo^p is the typical optical depth of a den- 
sity scale height. For r^m <C 1, the radiation has ample 
time to converge to the diffusive solution given by the in- 
stantaneous mass configuration. In the opposite limit, the 
radiation does not have enough time to diffuse and any 
perturbation to it is synchronized with the gas. 

Note that for the optically thin part of the atmosphere, 
one should consider the nonlocal nature of the diffusive 
solution. In such a case, the diffusive time scale is simply 
given by the time it takes the radiation to traverse a scale 
height, i.e., Tdiff ~ Ip/c, yielding rdiff « ct/c. 

A third time scale emerges from the energy equation 
(eq. 3). We look at the set of eqs. (l)-(5) in the linearized 
case, when zero order gradients are neglected and when the 
terms arising from the finite propagation speed of light can 
be neglected as well. After applying a Fourier transforma- 
tion, eqs. (4) & (5) become 

where Xv is the total extinction per unit volume and 
and El are the Fourier transforms of Si and Ei respec- 
tively. We find that for (P/JiXvKv) :$> I, Si - Ei ^ Si 
while for {k'^/3xvKv) < 1, 5i - « 5i(fcV3x^K^). If 
we now plug this result into the linearized and Fourier 
transformed energy equation (cq. 3) and use the first two 
hydrodynamic equations (eqs. 1 & 2) to simplify it, we 
obtain after some algebra (in a similar manner to §101 of 
Mihalas & Weibel Mihalas 1984) that 

{ [u^ - Cak^] uj + T-ii(fc) y - CTk^] } i) = 0, (24) 
where 



1 



12/3(7- 1)k„c 



1 + 3 



fc2 



(25) 



Clearly, when T^g],i ^ w, the wave equation obtained will 
be the isothermal wave equation while for r^^^j <JC lo, the 



adiabatic wave equation is obtained. Tcooi{k) is therefore, 
the typical time scale which takes a wave with a wavevec- 
tor k to cool radiatively. 

We can define a dimensionless ratio between the cooling 
time scale for a wave with a wavelength of the order of the 
scale height and dynamical time scale that is the time it 
takes the same wave to cross a scale height as: 



''cool — Tcoo\{k — ^Tvlp )/Tdyn 



12/3(7 - l)cK.vlp 



12/3(7-1) T^1^^7 XvKvll<S:l 
TfW(^^^o for XvKvl^^l- 



(26) 



(27) 



For rcool ^ Ij the cooling time scale is too long for the 
perturbations on scales of the scale height to have time 
to equilibrate to the radiation temperature. The gas is 
then close to the adiabatic limit. In the opposite limit of 
fcool ^ 1, the gas has time to reach thermal equilibrium 
with the radiation and if the radiation has no perturbation, 
it is close to the isothermal limit. Since the equilibrium of 
the radiation can in principle be non-isothermal (if there is 
a net flux in the system), this limit is more generally called 
the NAR limit (non-adiabatic reversible), since in this op- 
posite limit to adiabatic perturbations, the perturbations 
are again reversible. 

Fig. 1 depicts the typical value of the dimensionless pa- 
rameters in a typical Thomson atmosphere as a function of 
the optical depth from the top of it. Four distinct domains 
are seen in the Tcooi/rdyn vs. Tdiir/Tdyn plane. The vari- 
ous approximations are depicted in this plane. The figure 
illustrates the run of the physical conditions in a typical 
star as one moves from the inside outwards, namely from 
high to low optical depths. The effect of the ratio of ex- 
tinction to absorption on the path in the parameter plane 
is also shown. The discussion in this paper concentrates 
on the l.h.s. of the parameter plane, that it, on the re- 
gion in which the radiation does not evolve synchronously 
with the instantaneous gas configuration, i.e., when the 
radiation and gas cannot be considered as one fluid. 

4. FULL GLOBAL LINEAR ANALYSIS 

We begin with a complete global linear analysis. By 
changing the free parameters of the problem (F, boundary 
conditions, etc.) we expect to find the various instabili- 
ties that exist in this type of systems. Afterwards, we will 
analyze each one of them separately. 

4.1. Equations and Problem set-up 

The governing set of equations is given in §2. After 
linearization and Fourier transforming the horizontal di- 
rection and time we get the following equations for the 
transformed quantities: 



. P ., - d _ .15 

lU} = tkU + TT-W -|- W — T— po 

po oz po oz 



d . 

-tupow = --^P 

oz 



+Po 



9*p-i 



Po 



(28) 
(29) 



En 



P , P 

Xp — \-Xp — 

Po Po 



6 



N. J. Shaviv 




Fig. 1. — The connection between the ratio of the cooUng time scale to the dynamical time scale and the ratio of the diffusion time scale 
to the dynamical time scale for a few saniijle Thomson atmospheres, as a function of optical depth. The case depicted by a thick solid line 
corresponds to 7 = 5/3, c^/c = 10^'^, ft/x = 10^^ and T = 0.5 {0 = 1 in the optically thick limit). The two lighter solid lines correspond 
to changing the absorptive opacity from n/x = 10^'^ to 1 or to 10^*. The dashed line corresponds to changing F from 0.5 to 0.05. The four 
regions correspond to whether the radiation is in the diffusion limit or not (i.e., whether it behaves as two fluids or one coupled fluid), and 
whether the gas is in the adiabatic or NAR limits. Convection is efEective in the top half. The instabilities described here operate in the left 
half. 
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-iLoF, + j-Q^E = -po(k + (t)cF, - p{k + a)cFo (33) 
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where w is the vertical component of the velocity and u 
the lateral one. These cqtiations are the equations that ap- 
pear in (Spiegel & Tao 1999) with two differences. First, 
we write the radiation in the co-moving system with the 
material. Second, we retain the option of a variable opac- 
ity for generality and usage in the next work. We omit for 
brevity the '1' index from the Fourier transformations of 
our eight variables pi, pi, mi, wi, Si, F^.i and 
These variables satisfy four complex algebraic equations 
and four complex ODE's with the height z as the inde- 
pendent variable. 

The Boundary Conditions: The set of ODE's re- 
quires boundary conditions. Since the dependent variables 
p, w, E, Fz are complex variables, four complex boundary 



conditions arc required. 

Several gas boundary conditions can be imposed. At 
either the bottom or the top we choose m; = or p = 0, 
which result with reflection of the waves at the boundary. 
We find that these conditions do not change the qualita- 
tive behavior of the waves. We should note however, that 
if we wish to compare our results to the behavior of real 
atmospheres, the results will have any meaning only if the 
waves in the real atmosphere considered can be trapped 
(and not for example propagate outward where they can 
be dissipated). We shall elaborate this point more in the 
discussion section. 

As for the radiation, there are more possibilities. 
The following possibilities exist for the top: 

Fixed radiation temperature: E{zt) =0 or 

Eddington t = condition: E{zt) = 2F{zt)/c 
While the following possibilities exist for the bottom: 

Fixed radiation temperature: E{zb) = or 

Fixed incident flux: Fi^b) = 

When the top of an atmosphere is free to emit radia- 
tion to infinity, the Eddington condition for the top is the 
more appropriate condition. The fixed radiation tempera- 
ture will be the physically appropriate one if for example 
the layer of interest is bounded at the top by a layer that 
equilibrates quickly its temperature by convection. As to 
the bottom conditions, they too can depend on the system 
at hand. For example, if the atmosphere is above a layer 
which conducts heat well, its temperature will be fixed. 
In principle, this boundary condition can be a mix of the 
two. 
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4.2. Form of Solution 

In order to find the eigenmodes of the acoustic oscil- 
lations, wc first need to construct a profile of the unper- 
turbed atmosphere. We avoid at this stage the usage of 
real stellar atmospheres since we are more interested here 
in understanding the physical origin and behavior of the 
unstable modes in the simple Thomson scattering domi- 
nated case and not to quantify the unstable behavior of 
real stars that have complex atmospheric opacities. More 
realistic opacities as well as non LTE effects will be in- 
cluded in future work. As stated before, wc also neglect 
for simplicity the non-slab geometry of most relevant at- 
mospheres. 

The equations describing the unperturbed steady state 
atmosphere are the radiation diffusion equations: 



dEp 

dz 



const 



and the hydrostatic equation 



dpo 2dpo . , XvfiFo 



(36) 
(37) 

(38) 



Since the gas and radiation arc in thermal equilibrium in 
the unperturbed solution, wc also have = cE/AasB. 
which allows the calculation of c-r . 

Inasmuch as the Eddington relation should be satisfied 
in particular in the optically thin limit, we require that at 
the top of the atmosphere E{t ^ Q) = 2F/c. 

The main parameter of interest in this work is F = 
XvfiF/cg. Hence, we solve for an unperturbed atmosphere 
with a given F extending over a given optical depth range 
(say, between r = 0.1 — 100). The solution is carried out 
numerically and then tabulated so that a simple spline in- 
terpolation can be applied. The values of p, p oi E can be 
quickly and to any given prescribed accuracy be found for 
any height in the atmosphere. 

The next step is to find the eigenmodes of waves with 
a given horizontal wavevector kx in a given atmosphere. 
Our set of ODEs has four first order complex ODE's and 
complex eigenvalues. Two of the boundary conditions are 
given at the bottom of the atmosphere while the two oth- 
ers are given at the top, making it a complicated boundary 
value problem. To find an eigenmode we use the shooting 
method. The two 'missing' bottom boundary conditions 
are iterated for until the imposed top boundary conditions 
are satisfied. The procedure of finding the solution is car- 
ried out in the following way: 

1. Since two boimdary conditions are given at the bottom, 
we guess the values of the two other variables that don't 
have a specified value on the boundary along with a 
guess for the unknown eigenfrequency u. For example, 
if p = and E = are the bottom boundary con- 
ditions, we start by guessing the values of F^, w and 



2. Once all the dependent variables are given or guessed 
on the lower boundary, we can integrate upwards. The 
set of equations can be written as four ODE's and four 
algebraic equations that can be evaluated at each point 
z in the following order: 



3. 



Evaluate K,a,x = K + a,Eo,Sp = —po ^{dpo/dz) 
and g*=g - {x/c)F. 

Calculate numerically the variables with algebraic 
equations: 

kc'^E 



F^ = 



S = 



3{ucU} + ipoXc) 
kp 

U!po U) C 



J- X 



(39) 

(40) 



iujp - clpo{ico h Spw) + E{'y - 1)pockoWk ^ 

Po 

+pQg*w] I [(7 - 1)poc«;o/wk - iu^c^ IA.E^ (41) 



P 



S 



• Integrate upwards using the following ODE's: 

d . . p ., . . 

— w = tu tku + wsp 

OZ po 

d X 
-^p = tupow - g*p + po-Fz 



(42) 



(43) 



dz 
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+P0F0 
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Po Po 



(44) 



4:En 



-^Fz = UKpoiic{S — E) + UciioE — ikFx — iui— — Ud 



dz 



-3po-Fz + 



Sp^Fo 



-3po-Fo 
c 



P , P 

Xp — \-Xp — 

Po Po 



3po 
(45) 



(46) 



We marked several terms so as to be able to trace 
their effect. The terms multiplied by Uc are terms 
that originate from the wavy behavior of the radi- 
ation field, while terms multiplied by Uk are terms 
that originate from the finite cooling time of the sys- 
tem. The terms with are the 'Doppler' terms. In 
the normal case Wc = = Wd = 1- In the limit 
Mc ^ wc get the instantaneous limit for the radia- 
tion. While in the limit u^, we recover the negli- 
gibly small gas heat capacity limit (the NAR limit). 
The limit Ud ^ shuts off the effects of the Doppler 
correction in the radiation field to the motion of the 
material. Since one of the goals is to check whether 
these terms are important for instability, they are 
specifically marked and can be artificially changed, 
though they are equal to unity in real cases. 

The numerical integration is accomplished using the 
fourth order Runge-Kutta integration scheme. 

At the top, the solution is compared with the two top 
boundary conditions. The initial guess for the two ad- 
ditional bottom conditions and the eigenfrequency iv 
are then corrected and the equations are integrated up- 
wards once again. The 'shooting' upwards it iterated 
until the solution converges to the required top bound- 
ary conditions. The solution obtained is an eigenmode 
of the system and its eigenvalue is the converged fre- 
quency CO. 
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Since we are interested in finding more than one partic- 
ular eigenmode, we use the above procedure to search for 
more modes. This search is performed by systematically 
guessing different initial guesses for the eigenvalue w. In 
this way, the entire frequency domain is covered up to a 
given frequency. 

In principle, other types of modes such as radiation dif- 
fusion modes at very high frequency can be searched for 
and found, however, they are beyond the interest of this 
work as they arc generally very stable (with large damp- 
ing rates, e.g. Mihalas & Weibel Mihalas 1984). Sim- 
ilarly, we are not interested in convection or g-modes. 
For any reasonable absorptive opacity larger than ^ 
(ca/c) (prad/Pgas)c, the top part of the atmosphere will 
cool too quickly for convection to be efficient. This will 
appear in the form of modes that have are real (3-modes) 
or pure imaginary (convective modes) with absolute os- 
cillation or growth rates that are much smaller than the 
dynamical time scale. Only if we force the atmosphere 
to be adiabatic by artificially having a very small absorp- 
tive opacity, will the jf-modes or convective modes obtain 
a 'dynamic' time scale. 

Once the eigenmodcs of waves with a given and F arc 
found, we can systematically search for the eigenmodes of 
other kx, or of atmospheres with a different F. Since we 
are interested in the locus of the eigenvalues as a function 
of either kx or F, we can evaluate the locus by changing 
kx or F by small increments and use the previous value of 
the two independent variables and eigenvalue as the new 
initial guess. We assume that the eigenvalues are contin- 
uous functions of k^ and F. This allows a relatively fast 
construction of the locus of eigenvalues. 

4.3. Results 

We divide the description of the results into three parts. 
We first begin by analyzing the eigenvalues u) of an at- 
mosphere with a given F and horizontal wavenumber kx 
and summarize the two instabilities found. We proceed to 
construct the eigenvalue spectrum of an atmosphere as a 
function of the horizontal wavevector kx and then perform 
a similar analysis to find the eigenvalues as a function of 
the Eddington parameter F. In subsequent sections, we 
will then study the two absolute instabilities. 

4.3.1. An atmosphere with a given F and kx and the 
instabilities found 

As described in §4.2, for a given set of boundary con- 
ditions imposed, the eigenvalues corresponding to a given 
F, kx pair are found. 

Fig. 2 depicts the lowest acoustic eigenvalues obtained 
for three atmospheres with F = 0.9 and kx — ^/Ip- One at- 
mosphere has a fixed temperature imposed at the bottom 
boundary condition while the other two have a fixed flux 
imposed at the bottom. In all three cases the Eddington 
approximation condition is imposed at the top (namely, 
Et = 2Ft/c), at Tt = 0.3. In the fixed temperature case, 
the bottom is fixed at = 10. In the two fixed flux at- 
mospheres, the bottom is located at Tf, = 10 and 30. 

We assume that there is a small absorptive opacity 
K = x/100. This absorptive opacity, however small, is still 
100 times larger than the value needed to force the per- 
turbations with a period of the dynamical time scale to 



be in the NAR limit and not the adiabatic one. This is a 
better approximation for stellar atmospheres that most of- 
ten have a small but significant absorptive opacity in this 
respect. One the other hand, the absorptive opacity is 
too small to change the total opacity which is dominated 
by electron scattering. The gravitational acceleration in 
this particular case is lO^cm s~^, an acceleration which 
corresponds to a White Dwarf envelope. The physics ob- 
tained is to a large extent not a function of this parameter. 
Choosing a different gravitational constant will just imply 
that the typical time scale will be different and the effec- 
tive temperature needed to sustain a given F is different 
as well. In this particular case, the effective temperature 
obtained is 5.82 x 10^ °K. At this temperature, the ra- 
tio ct/c at the center of the layer is roughly 1.0 x 10~^. 
(In reality these conditions are similar to those found in 
novae.) 

We can now sec the two main instabilities. When the 
lower boundary condition has a fixed temperature, the low- 
est vertical mode can be unstable to the first type of in- 
stability - 'Type F. When unstable, this mode's eigenfre- 
quency becomes purely imaginary. It grows on dynamical 
time scales while its counterpart decays with the same rate. 
The growth rate found in this case is 7.7s~^. This is the 
dynamical time scale of the system. 

Under the two different lower boundary conditions, an- 
other type of instability can operate - 'Type IF. When a 
mode succumbs to this instability, is has a real part with a 
somewhat smaller but significant positive imaginary part. 
Namely, it is mostly traveling but it amplifies (and its 
counterpart decays) on a dynamical time scale. In the case 
depicted here w = (8 -I- 0.45i) s^^ with only small changes 
between the three cases. Although the growth rate of this 
mode is not as fast as fast as that of the Type I instability, 
the e-folding growth time of this mode is nevertheless only 
three oscillation periods. Moreover, unlike the Type I in- 
stability, this one is relatively insensitive to the boundary 
conditions. 

Fig. 3 shows the structure of both unstable and stable 
modes. The main features seen are that the Type I mode 
has gas and radiation energy densities that are inversely 
correlated with each other. In the Type II mode, there is 
clearly a phase lag between the density and the radiation. 
This lag flips sign in the conjugate pair of the Type II 
unstable mode. 

We also checked whether the various artificial factors 
that we introduced change the occurrence of the modes. 
Uc changes the ratio between the speed of light and speed 
of sound. We found that both instabilities exist in the 
two fluid limit (in which the radiation speed is infinite) 
but disappear when in the one fiuid limit, when the ra- 
diation cannot diffuse. u„ changes the importance of the 
absorption opacity which determines the ratio between the 
cooling time scale and the dynamical time scale. Wc found 
that both instabilities exist in both limits, the adiabatic 
limit in which the gas has no time to equilibrate with the 
radiation temperature and in the NAR limit, in which the 
gas temperature is set by the radiation temperature. 

We also searched for a third instability like the one de- 
scribed by Spiegel & Tao (1999) but could not find it. We 
defer further discussion on it to the appendix where we 
will try to understand how such an instability could arise 
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artificially. We also explain why the two dynamic instabil- 
ities found here were not found by Spiegel and Tao even 
though they did solve the same equations. 

4.3.2. Eigenmodes as a function of kx 

Taking the two examples depicted in Fig. 2 with = 10, 
we now look at the locus of eigenmodes obtained as a func- 
tion of kx- This is depicted in Fig. 4 for the fixed temper- 
ature case, which has both Type I and Type II unstable 
modes, and in Fig. 5 for the fixed flux case, which has only 
Type II unstable modes. 

The Type I modes are found to be unstable for all kx 
smaller than a critical wavenumber. This wavenumber de- 
creases as r is reduced. That is to say, once the critical 
Tcrit = 1/2 for instability is surpassed, it is the longest 
wavelengths which first become unstable. However, since 
3(a;) is found to be proportional to kx at long wavelengths, 
the most unstable mode has a finite horizontal wavelength. 
This is seen in the additional lines in Fig. 4 which are la- 
beled with varying F's. The additional lines describe the 
growth rate of the Type I unstable mode for different F's 
with the most unstable horizontal wavelength marked with 
a small open circle. 

The Type II modes are seen to occur only in finite wave- 
lengths. In Fig. 5, we can actually see the second eigen- 
mode become Type II unstable in one region of kx by 
merging with the third eigenmode and in another region, 
by merging with the first eigenmode. As the F is reduced 
towards the critical F needed for instability, we see that 
the most unstable mode does not vary by much. The in- 
stability is triggered at a finite horizontal wavelength. It is 
found to be triggered at Fcrit ~ 0.86 and a wavenumber of 
fcx.crit ~ 0.88 (?~^). These numbers are found to vary by 
a few percent if^ either the absorption opacity or the depth 
of the atmosphere are changed. 

4.3.3. Eigenm,odes as a function o/F 

Last, we study the behavior of an atmosphere as a func- 
tion of the Eddington parameter F. We take the example 
depicted in Fig. 4 and instead of changing k^, we keep it 
fixed at kx = (ip^) and change the Eddington parameter 
from F = 0.1 to F = 0.98. 

It is clear from Fig. 6 that Thomson scattering domi- 
nated atmospheres become unstable at large enough lumi- 
nosities. However, the critical luminosity is well below the 
Eddington limit. They range from 0.5 in the fixed tem- 
perature case to about 0.86 when a fixed flux is imposed 
at the bottom. 

The main feature apparent in the flgure is that as the 
radiative flux increases, more and more modes become un- 
stable to the Type II instability. That is to say, it isn't 
only the lowest or second lowest vertical eigenmodes which 
are unstable. The dynamic state that this atmosphere will 
reach obviously depends on the nonlinear behavior of the 
modes and their mutual interaction. It is clear however 
that it will not stay homogeneous close to the Eddington 
limit. 

5. THE TYPE I INSTABILITY 

The two main clues pointing to the origin of the Type 
I instability were seen in the previous sections. First, is 
was found that the instability arises only when the bottom 



boundary condition has a fixed temperature, as opposed 
to a fixed fiux. Second, when it does arise, there is clearly 
an anti-correlation between the radiation energy density 
and the material density. We now proceed to show how 
this instability can arise with a toy model. 

We begin by showing how the boundary conditions af- 
fect the correlation between gas and radiation. We will 
then proceed to show how an anti-correlation results with 
the Type I instability. 

5.1. Non local effects arising from a global analysis 

Since the boundary conditions were found to be cru- 
cial, we evidently cannot study the instability with a local 
analysis. We begin with an atmosphere in the NAR limit. 
Thus, at any given instant, the radiation field is governed 
by the equilibrium diffusion solution: 

V-F = and F = ^^, (47) 

and the boundary conditions imposed at the bottom and 
top of the atmosphere. 

To demonstrate how a positive or negative correlation 
can arise between Ei and pi, we explore a specific ex- 
ample. Suppose we have an unperturbed atmosphere in 
which xo — const that is perturbed by 5xv = Sxv{x), i.e., 
that the perturbations are only a function of the horizontal 
direction x. We also assume for simplicity that the verti- 
cal variations in the atmosphere can be neglected at this 
point. This implies that the flux F and its perturbation 
do not depend on the vertical coordinate. 

The optical depth between the top and any given point 
in the atmosphere is 

T=/" Xvdz = To{l + 5xvix)/Xv,o) , (48) 

J z 

where to{z) is the unperturbed optical depth from infinity 
down to the mass element at z. 

Next, we need to impose boundary conditions. If we 
impose r = at the top of the atmosphere, then from 
the Eddington approximation we have that at the top 
Et = 2Ft I c. At the bottom we can impose a fixed temper- 
ature and therefore a fixed i?;,, or, we can impose a fixed 
flux i^fc. An example for the implementation of a flxed 
temperature (or flxed E}f) at the bottom is the case of a 
very good energy conducting layer below the layer under 
consideration. 

The fixed flux boundary condition at the bottom cor- 
responds to the case when the layer below the one under 
consideration has a very large heat transfer resistivity so 
that the gradient in the layer sets a fixed flux that cannot 
be changed by the atmosphere exhibiting the instabilities. 

When the temperature and energy density E\, are cho- 
sen as fixed on the lower boundary, the flux F and the 
energy density at any height are found to be; 
cE^, 



F 



E = 



2 + Zn 
2-H3T 



= const 



Eh 



2 + in, 

Eo{to) (l - 



G(7 



2 + ST{l + Sxvix)/xv,o) 

2 + STb 



(49) 



(2 + 3t6)(2 + 3to) 

where Eo{tq) is the radiation energy density at the point 
where the unperturbed optical depth was tq. This result 
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Fig. 2. — The eigen frequencies of the lowest modes in three atmospheres. The two panels refer to the two possibilities of the sign of rj 
— the growth or damping rate of the modes. The filled circles are the modes of an atmosphere with F = 0.9, kx = i-/lp, g = lO^cm s~^, 
k/x = 1/100 and a fixed temperature imposed at a bottom placed at rj, = 10 while the top is at r = 0.3. The empty circles are the modes of 
the same atmosphere with a fixed flux imposed at the bottom. The crosses are the eigenmodes obtained when the bottom boundary of the 
second atmosphere is moved to t = 30. The Type I instability is seen in only the first atmosphere with a fixed temperature at the bottom. 
The Type II instability is seen in all cases. 
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Fig. 3. — The profiles of unstable and stable modes. The top left panel describes the Type II unstable mode of the second atmosphere of 
the previous figure (fixed flux at the bottom). The bottom left panel depicts the profile of the conjugate mode which is highly damped. The 
top right panel shows the profile of the Type I unstable mode of the first atmosphere described in the previous figure (fixed T at the bottom). 
The bottom right panel shows the profile of the fourth mode in the same atmosphere. It is slightly radiation damped. The variables plotted 
are the energy density E (dashed lines) and the matter density p (solid lines) multiplied by /3cy, so as to be on the same par as E. Since the 
problem is linear, the absolute normalization of the modes is meaningless. The real part is given by heavy line while the imaginary part by a 
thin one. 



implies that under a positive perturbation to the opacity, negative: 

the change 6E{tq) in the radiation energy density wiU be 5E{to) 6(rfc — tq) Sxvix) 

E ^~(2 + 3Tfc)(2 + 3To) ' 



(50) 
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Fig. 4. — The eigenvalues of the first atmosphere depicted in Fig. 2 that has a fixed temperature imposed at its bottom, as a function of the 
horizontal wavenumber kx in units of the pressure scale height. The bottom panel is the real part of the eigenmode while the top panel is the 
logxo °f absolute of the imaginary part. Unstable modes are denoted by thick lines. This atmosphere exhibits both Type I and Type II 
unstable modes. Below the critical horizontal wavenumber, the Type I unstable mode has a purely imaginary eigenvalue. Above the critical 
horizontal wavenumber, the absolute value of the real part is much larger than the absolute value of the imaginary part namely, there is a 
small damping term arising from radiation diffusion. The additional double lines represent the imaginary part of the Type I instability when 
r is decreased from 0.9 to 0.7, 0.6, 0.55 and 0.51, with which we see that the most unstable wavenumber moves to the long wavelength limit. 




Fig. 5. — Same as the previous figure except that now the flux is kept constant as the bottom boundary condition. Here, the atmosphere 
is only unstable to the Type II instability. We see that the second eigen mode 'merges' at two different horizontal wavelengths, with either 
the third eigenmode or the first eigenmode. The additional lines represent the imaginary part of the growth rate as F is decreased. We see 
that the most unstable wavenumber changes only slightly with F. The growth rate is typically several times the oscillation period. 



that is to say, there wiU be an anti-correlation between 
6E and Sxv The coefficient of the anti-correlation isn't 
large. The largest anti-correlation is obtained if the lower 
boundary condition is fixed to be at Tb — 2\/2/3 ~ 0.94 in 



which case SE/E = -{'S-2V2)Sxv/Xv,o ~ 0.172^Xf /Xf ,o- 
The anti-correlation is larger than 0.1 for 5 > Tf, > 0.2. 

A negative correlation implies that the radiation pres- 
sure is anti-correlated with the density. 
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Fig. 6. — The eigenvalues of the first atmosphere described in Fig. 2 (with a fixed T at the bottom) as a function of F, while kx is kept fixed 

at (ip^^ . uj Ss multiplied by 1/(1 — F) to normalize it to the changing dynamical time scale as the atmosphere puffs up when F is increased. 

The graph is actually plotted as a function of logiQ(F/(l — F)) « logjQ P{t » 1) in order to extend the high F region. Like in the previous 
figures, unstable modes are depicted with heavy lines. We see that as the flux is increased, this particular atmosphere first becomes unstable 
to the Type I instability (with a purely imaginary mode). As F approaches the Eddington limit, more and more modes become unstable to 
the Type II instability (which are mostly propagating). 



When a fixed fliix is imposed at the bottom, we have 



F = Fq = const 
E = {2 + 3t)^ ^ Eo{ro) { 1 



3to Sxv{x) 
2 + 3to Xv,o 
Here the correlation is positive, namely: 
6E{to) 3to Sxv{x) 



E 



2 + 3ro Xvfi 



.(51) 



(52) 



The coefficient of the correlation can be in this case larger 
than in the previous one and it increases with the opti- 
cal thickness of the atmosphere. For Tfe ^ oo, we have 

6E/E = 6xAx)/Xvfi- 

We conchidc that under the same local conditions one 
can obtain qualitatively different behaviors according to 
the boundary conditions imposed. Should we have eval- 
uated this atmosphere using a local analysis, wc would 
have obtained that (^Prad) = since kz = 0. Clearly, 
the boundary conditions are important to the question of 
stability. 

In a real atmosphere with vertically dependent unper- 
turbed state and modes, wc could expect to sec cigcnmodcs 
with both types of correlations. We point that the correla- 
tions do not depend on /3 and hence such phenomena can 
occur in low luminosity atmospheres. A typical example is 
the formation of grains in atmospheres. The appearance 
of grains causes large opacity changes. 

5.2. Toy model with non local effects 

Wc assume for simplicity that the relevant perturbations 
(of the order of the scale height of the atmosphere) are in 
the NAR limit. We also assume that the gas is held within 



a narrow slab without any z dependence. That is to say, 
wc assume that the perturbation to the radiation density 
is proportional to the density perturbations: 

S p S p 

El = sEq — or equivalently Prad,i = ePrad.o — , (53) 
P P 

where e is the 'efficiency' with which the average radia- 
tion energy density or pressure changes with changes in 
the opacity (which is a function of p). In §5.1, e was cal- 
culated for a very simplified case. In reality, ^Prad is not 
simply proportional to Sp but is also a function of height, 
we just assume it to be so (this is why it is a 'Toy' model!). 
We therefore leave s as an effective free parameter. 

The equations describing the gas are the continuity and 
momentum equations. The linearized versions of which 
are 

dp._ ...... (54) 



upi „ , 
— = -Vpo • vi - (V • vi)po 



9^1 „ „ 

PO-^ = - Vpi - Vprad,: 

at 



(55) 



Since the gas is in the NAR limit, its pressure perturbation 
is given by 

El 



2 

Crppl. 



(56) 



The radiation pressure perturbation can be written as 

Prad,l = PPO— = Pc^Pl- (57) 

Po 

Thus, 

^ A/91 (l + £/4 + 5/3)4 = 0. (58) 



9*2 



Radiative Instabilities 



13 



The first term in the parenthesis arises from the isothermal 
speed of sound and describes simple acoustic waves. The 
second term arises because the NAR limit is not necessar- 
ily the isothermal limit. Consequently, the residual corre- 
lation between 5T and 5p imply that the relevant sound 
speed is modified to: 

4AR = 4(l + e/4). (59) 

The third term within the parenthesis arises from the work 
that the radiation does in synchronization with the wave. 
The corresponding dispersion relation is 

= fc2 (1 + e(/3 + 1/4)) 4 = k'^cl^. (60) 

A negative e implies a decrease in the effective speed of 
sound. When the anti-correlation is large enough, the 
speed of sound c^g- becomes imaginary and the waves sim- 
ply amplify without propagating. In most astrophysical 
systems, this cannot arise unless /3 is sufhciently large be- 
cause e is seldom large by itself^. 

If e is negative, the system will necessarily become un- 
stable at some L ~ ^crit < ^Edd for which (3 is large 
enough. Using the toy model and a typical value below 
the photosphere of £ « —0.15 which we should expect to 
see in an atmosphere with a Tf, of about 5, we have that the 
critical /? in this case is /3ciit ~ 1/e— 1/4 « 6.5, and the crit- 
ical Eddington parameter is Fcrit ~ /3crit/ (/^crit + 1) ~ 0.86. 
The full numerical calculation gave Fcrit = 1/2 for a deeper 
atmosphere. This discrepancy can actually be reconciled 
exactly with a more elaborate toy model. The feature 
missing in this simple model is the 'feedback' that the 
change in the radiation flux has on the vertical extent of 
the atmosphere. Specifically, a lower density region will 
have a higher flux going through it and therefore a smaller 
effective gravity. Consequently, it would puff up. At a 
given height above the bottom, the radiation density will 
therefore be even higher because the Lagrangian coordi- 
nate will be lower. This will introduce a further anti- 
correlation between p and E. This missing feature in the 
toy model also allows very deep atmospheres to be unsta- 
ble. Without it, we saw in §5.1 that as Tb is increased, the 
anti-correlation tends to zero. 

Once an atmosphere is unstable to Type I instabilities, 
it would tend to open 'chimneys' through which it is easier 
for the radiation to escape and also to accelerate material. 

6. THE TYPE II INSTABILITY 

To understand the Type II instability, it is best to 
compare it to the instability of radial s-modes (e.g., 
Glatzel 1994, and Papaloizou et al. 1997). As we shall 
see, there are many similarities and a few particular dif- 
ferences which will elucidate the origin of this instability. 

The main similarities between s-modes and the Type I 
instability are: 

1. The instabilities result with eigen- frequencies that de- 
scribe modes with e-folding growth times of order the 
oscillation period. Each mode has a conjugate that is 
damped at the same rate. Moreover, as a control pa- 
rameter is changed (e.g., Tes in s-modes or F here), 
the conjugate pair arises when two modes with differ- 
ent real frequencies merge together. 

^ However, the following axe two known exceptions, (i) A photon + 
with grains in which the opacity changes can be extremely large. 



2. Both types of instabilities appear only in systems in 
which the radiation pressure dominates. 

3. Both instabilities are quite different from various 'slow 
instabilities' such as the K-mechanism and other Carnot 
types. 

4. Both instabilities disappear when the system is driven 
out of the two fluid limit. Namely, both need the ra- 
diation field to be described by the diffusion equations 
instead of being highly coupled to the gas. 

Nevertheless, the are a few critical differences between 
the two: 

1. The Type II instability operates in Thomson atmo- 
spheres, s-modes need a particular opacity law for them 
to be unstable. In particular, they need kt > which 
is obtained in ionization zones. Thomson atmospheres 
are much more general. 

2. Unstable s-modes are a radial phenomenon. The Type 
II is intrinsically non-radial. Type II unstable modes 
do not exist if the horizontal wavelength is significantly 
larger than the pressure scale height of the system. It 
also implies that it will be a local phenomenon in stars 
since the scale height is usually much smaller than the 
radius. 

3. s-modes are localized to the top part of the atmosphere, 
where the NAR limit is obtained, by having a region 
with an effective speed < 0. The Type II unstable 
modes will be localized to the top because they are a 
very high i phenomenon. 

4. Last, unstable s-modes require that the gas tempera- 
ture be given by the radiation temperature, namely, the 
NAR limit. Type II unstable modes where found to ex- 
ist also when the atmosphere was switched to the adia- 
batic limit (though still in the two fiuid limit) by reduc- 
ing the absorptive opacity enough to have TcooI ^ Tdyn- 

We sec that, mathematically, the origin of the two dif- 
ferent instabilities is similar. The instabilities arise when 
the linearized equations become non self-adjoint (e.g., Pa- 
paloizou et al. 1997) as it is this characteristic that results 
in the complex conjugate pair of eigenmodes. In the case 
of s-modes, it was shown that it is because the dispersion 
equation becomes third instead of second order. Without 
the interaction with the diffusive radiation field, the pres- 
sure and density have a local algebraic relation of the form 
Pi = v^pi. As such, the dispersion equation for is self- 
adjoint and no complex roots are found for ui^. However, 
because the temperature is set by the radiation field, the 
relation between density perturbation and total pressure 
perturbation becomes differential (Glatzel 1994). To see 
this, we look at the vertical component of eq. (5), and 
perturb it: 

dz dz 3 dz c c 

gas fluid with i^'s playing the role of the radiation, (ii) An atmosphere 
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In the case of radial perturbations, F^^i vanishes. Thus, 
in s-modes, it is the first term on the r.h.s from which 
a differential relation is obtained between ptot,i and pi 
which Xv,i is a function of. The dispersion relation be- 
comes third order and cj^ can obtain complex conjugate 
roots. The differential nature of the pressure density re- 
lation is obtained by perturbing the radiation transport 
equation under a fixed flux, as it is in the two fluid limit 
(Glatzel 1994). Because the variable opacity, which is a 
function of the gas temperature, is one of the key ingredi- 
ents needed to get the differential relation, it is clear why 
s-modes do not exist in the two fluid adiabatic limit. In 
this limit, the temperature of the gas, and therefore the 
opacity, is not a function of the radiation field but of the 
gas alone. 

The physical origin of the Type II instability is similar 
to that of s-modes though not the same. The diff'erent 
characteristics of the two instabilities should point to the 
differences in physical origin. 

First, like s-modes, the Type II doesn't exist in the 
one fluid limit and the radiation pressure must dominate, 
thus, it is probably the diffusive nature of the radiation 
field which promotes the Type II instability. Second, since 
Type II exists also in Thomson atmospheres, the form of 
the opacity is obviously not crucial. It also explains why 
the Type II exists also in the two fluid adiabatic limit - 
the instability is not sensitive to the properties of the gas, 
whether its tcmpciraturc is set adiabatically or by the ra- 
diation. Thus, we hypothesize that the origin of the Type 
II instability is in the differential nature that the relation 
between the radiation pressure perturbation and the den- 
sity perturbations obtains due to non-radial perturbations. 
The second term in cq. (61) fits the required features. It 
will not vanish if F^^i ^ 0, namely, if there are perturba- 
tions to the flux in the vertical direction. Since V-F = in 
the diffusive limit, such a perturbation to F can arise only 
if there are horizontal perturbations F^.i as well. These 
arise only when there are horizontal perturbations to the 
density on the same scale as the vertical perturbations, 
hence the need for non-radial modes. 

To check this hypothesis, we artificially modify in the 
radiation equations the term responsible for horizontal 
fluxes. Namely, we write VE {dE/dz)z + Uh{dE/dx)x, 
where is an artificial factor. For u/j = 1 we recover the 
correct diffusion equations. For Uh ^ we inhibit the flow 
of radiation in the horizontal direction. If our hypothesis 
is correct, then under the latter limit, the Type II insta- 
bilities should disappear. What we find in the numerical 
calculation is just that. 

To summarize, the Type II instabilities originate from 
a differential relation between the total pressure and the 
gas density. This transforms the dispersion equations to 
a higher order equation in u) that can have complex con- 
jugate roots. The differential 'equation of state' needed 
to relate the density and total pressure perturbations in s- 
modes, arises when perturbing the constant flux equations 
and it is sensitive to the opacity. The required differential 
relation in Type II unstable modes arises when horizontal 
perturbations are taken into account, allowing horizontal 
fluxes. 

How does a Type II unstable system look like? To see 
this, we plot in Fig. 7 a snapshot of the eigenmode of the 



system. Plotted is the energy density of the radiation as 
a function of x or — t. The solution is proportional to 
Ei{z) exp{ikxX — i^{co)t) exp(9(a;)f). Thus, if we fix t, we 
obtain Ei which varies with z and is harmonic in the x 
direction. If we vary the snapshot with time, we will find 
that it propagates to the right side. If we fix x, let t and 
z vary, and exclude the exponential growth, we will ob- 
tain the same figure, with the horizontal axis for the time 
flipped. Namely, we will see at a flxed x, structure mov- 
ing downwards. This is opposite to conventional photon 
bubble picture. 

7. DISCUSSION 

We have seen that atmospheres with a significant ra- 
diative support always become unstable if the radiative 
pressure is increased sufficiently and becomes dominant. 
Both instabilities found are intrinsically non radial. Ta- 
ble 1 summarizes some of the characteristics of the two 
instabilities found and emphasizes the differences between 
them and other known instabilities - convection and s- 
modes which are both dynamic instabilities, and the k- 
mechanism which is a Carnot type instability. 

Using a Toy model, we can understand the Type I insta- 
bility. The anti-correlation between gas density and total 
pressure, dominated by the radiation, reduces the effective 
speed of sound. If the reduction is large enough, the speed 
of sound squared becomes negative and we have a purely 
imaginary frequency. This explains why 3(u;) was found 
to be proportional to kx at long wavelengths. 

We have also seen with the Toy model that in addition 
to the radiative pressure term, changes to the tempera- 
ture can contribute towards an instability. In Thomson 
atmospheres, this term is never large enough. However, 
when opacity variations can be very large, as is the case 
with dust formation, the temperature perturbation can by 
itself make the speed of sound become imaginary. In such 
a case, the lower limit for an instability is the lowest flux 
for which oscillations on the dynamical time scale are still 
not adiabatic. 

The 'temperature' term contribution to the effective 
speed of sound can become important in a second type 
of cases in which the equation of state becomes signifi- 
cantly 'softer' than that of an ideal gas. That is, when 
d\np/dlnT\p » dlnp/dln p\t- For example, in a system 
in which the fluid is composed of highly coupled radiation 
and plasma and the role of the radiation is played by a 
strong u flux, then given temperature variations will now 
be more important in setting the effective speed of sound. 
When this happens, we can again expect to see an insta- 
bility even in systems which are far from the Eddington 
flux (as long as oscillations on dynamical time scales are 
not adiabatic). 

The second type of instabilities found, the Type II, were 
found to be similar to s-modes. The difference however 
is that Type II does not require special opacity laws on 
one hand, but on the other hand, it does require non ra- 
dial oscillations. The latter allows for a differential re- 
lation between the radiation pressure perturbations and 
the gas density which can result with a non self-adjoint 
dispersion equation for w^. Once triggered, the unstable 
modes describe highly distorted horizontally propagating 
waves as is depicted in Fig. 7. Interestingly, if a standing 
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k^x/2TT or -Rc(co)t/2TT 

Fig. 7. — The radiation energy density of the Type II unstable mode of the 2""^ atmosphere of Fig. 2 as a function of z and either x or 
—t. If we choose the horizontal axis to be x, we see a snapshot of the eigenmode. Solid lines are positive iso-contours of the radiation energy 
density while dotted are negative. The normalization is meaningless in the linear problem. With time, this pattern moves to the right and 
amplifies according to exp(r]t). If we look at the evolution of the profile along a vertical cut, its evolution of E{z,t) without the exponential 
growth will be given by the same plot with the horizontal axis being —t, namely, it would appear as if structure is moving downwards. 



Table 1 

Comparison between the instabilities found here and three others. 



Instability: 


Type I 


Type II 


Convection 


s-mode 


K- mechanism 


Radial (exists for kx 0) 








+ 


+ 


Propagating / 0): 




+ 




+ 


+ 


Dynamic time scales (^(0)) ~ ©(l/rdyn)) 


+ 


+ 


+ 


+ 




Exists in adiabatic limit (rcooi/Tdyn ^ 00) 


+ 


+ 


+ 




a 


Exists in NAR limit (rcooi/Tdyn 0) 


+ 


+ 




+ 


a 


Exists in diflusion limit (rdiff/rdyn ^ 0) 


+ 


+ 


+ 


+ 


a 


Exists in one fluid limit (rdiff/rdyn 00) 






+ 




a 


Requires large prad/pgas 


+ 


+ 




+ 




Requires special opacity laws 








+ 


+ 


Significantly affected by boundaries 


+ 











The K-mechanism per se requires that the gas will not be fully in the adiabatic limit. Other similar Carnot type instabilities can 
exist when the gas is not fully in the other three limits. Namely, Carnot type instabilities can potentially exist if the gas is not in 
one of the four corners of fig. 1. 



wave in the horizontal direction is composed out of two 
oppositely propagating waves, the result appears as 'bub- 
bles' propagating downwards! That is, it appears as 'anti- 
bubbles'. The common conception that a 'bubbly' phe- 
nomenon should arise as the Eddington limit is approached 
assumes that the 'vacated' bubble through which it is eas- 
ier for the radiation to escape can be held together by 
something (Spiegel 1977). This is true in the case of 'pho- 
ton bubbles' in highly magnetized media where the ma- 
terial is forced to move only radially by vertical magnetic 
fields (Arons 1992). When no magnetic field is present, 
nothing is found to hold the material (Spiegel 1977), and 
as a result, these bubbles probably do not exit at all. What 
we find instead is a very dynamic picture in which the 



structure or 'phase speed' is actually downwards, such that 
the material in the 'vacated' regions is constantly changing 
on dynamical time scales. This is the opposite of bubbles. 
Of course, this pattern movement downwards does not in- 
volve net motion of material, and it facilitates the transfer 
of radiation upwards. What the exact form that this 'anti- 
bubbles' will look like, when in the nonlinear regime, is of 
course an interesting open question that requires a full 
numerical simulation. 

A very important question which should be addressed in 
future work is the effects that the assumption of trapping 
of modes has. That is, it was assumed in the form of the 
boundary conditions, that acoustic waves are reflected at 
the boundaries above and below the layer. This need not 
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necessarily be the case. If the atmosphere does not have 
a corona for example, where the temperature and speed 
of sound increase with height, any sound wave traveling 
through the atmosphere will just continue upwards and 
dissipate itself after forming shocks (which could in prin- 
ciple generate the required corona). In any case, we expect 
both instabilities to be important because they amplify on 
a short dynamical time scale - a mode does not need to 
travel back and forth to be amplified significantly. More- 
over, the imaginary effective speed of sound in the Type 
I instability will be shown in a subsequent publication to 
be a signature of a phase transition in which the homoge- 
neous state of the atmosphere is no longer the preferred 
or lowest energy state, that is, the instability found is just 
a manifestation of the system trying to get into a new 
equilibrium condition which is not the homogeneous one. 

Although a lengthy description of the nonlinear behav- 
ior is part of a forthcoming publication, a few words are 
in order about how such systems are expected to behave 
once these instabilities start to set in. 

The results of §4.3.3 have shown that either the Type I 
or the Type II instabilities will necessarily set in before the 
Eddington limit is reached. Thus, as a luminous object ap- 
proaches the Eddington limit, convection will become im- 
portant in the inner parts of the object where the density is 
high enough and the cooling time scale long (Joss, Salpeter 
& Ostriker 1973), while the outer parts will be unstable to 
the Type I or II instabilities where the two fluid limit is 
relevant, even if the opacity is Thomson dominated. Once 
either one of the instabilities takes place, it is clear that the 
atmosphere will not remain homogeneous. If the nonlinear 
pattern that is eventually formed propagates, as is the case 
in the Type II instability, each gas element will experience 
periodic density variations. While the peak force experi- 
enced by each mass element may exceed the Eddington 
limit, the time average force may be significantly below it 
(Shaviv 1998). In this case, an atmosphere can transfer 
through it a luminosity greater than the Eddington lumi- 
nosity without blowing itself apart. This will occur in the 
region of the atmosphere where perturbations of order the 
scale height become optically thin. From this region up- 
wards the effective opacity is again the microscopic one 
and the luminosity becomes effectively super Eddington. 
A predictable wind will then be blown from that region 
upwards (Shaviv 2000b). 

In the case of the Type I instability, on the other hand, 
the pattern formed is stationary and it resembles that of 
'chimneys'. The gas elements present in the lower density 
regions can witness a constant (in time) super-Eddington 
fl-ux. Although the average over the entire atmosphere can 
be a sub-Eddington luminosity, some gas elements in this 
case can be accelerated upwards through the chimneys, 
driving turbulence in the shear layers along the sides of 
the chimneys. However, this material will probably not 
be able to escape to infinity because it is unlikely that 
the escape speed will be attainable in the chimneys. This 
will require strong shears that presumable dissipation by 
shocks will prevent it from occuring. Thus, once the Ed- 
dington flux is surpassed in the chimneys, it will likely 
generate 'fountains' in which material accelerates upwards 
until the top of the chimneys is reached, then the material 
will fall back. If the average luminosity is super Edding- 



ton, the material need not stagnate and it can be easily 
driven off to infinity as a wind. 

In both cases, the clear role that the classical Eddington 
limit plays in the homogeneous case is lost. On one hand, 
a strong wind can exist already below the Eddington limit, 
driven for example b the energy dissipated from the un- 
stable modes; while on the other hand, quasi-stationary 
configurations (with a wind) can exist above the Edding- 
ton limit. 

A more detailed treatment of the nonlinear behavior, 
which includes toy models as well as numeral simulations 
is underway in a forthcoming publication. 

8. SUMMARY 

The two main results found in this paper are the exis- 
tence of two types of dynamically important instabilities 
that can take place in Thomson atmospheres. The two in- 
stabilities exist in addition to the possibility of convection. 
Following is a STimmary of the main results pertaining to 
these two important instabilities which will govern the be- 
havior of atmospheres shining close to the Eddington limit. 

1. Under a linear analysis of a slab Thomson atmosphere, 
two dynamically important instabilities were found to 

exist for the first time. 

2. Both instabilities operate in the two fluid limit of opti- 
cally thick atmospheres - when the diffusive time scale 
is shorter than the dynamical time scale of the sys- 
tem. Hence, they operate from the photosphere down 
to where convection can become efficient. 

3. The Type I instability is discovered once a nonlocal 
analysis of the radiation field behavior is carried out. 
A fixed flux boundary condition quenches it off. 

4. The Type I instability does not propagate. It corre- 
sponds to the opening of 'chimneys'. It appears in 
modes that have the radiation density on average anti- 
correlated with the gas density. The anti-correlation re- 
duces the square of the effective speed of soimd. When 
Cgg is negative, the instability appears and grows over 
a dynamical time scale. 

5. The Type I instability can take place also for 0, 
however, the growth rate then tends to vanish: r] oc 

0. 

6. A second instability, the Type II, appears under more 
general boundary conditions, but at somewhat higher 
F's than the Type I. 

7. The mathematical origin of the Type II instability is 
similar to that of s-modes. Since the hydrodynamic 

equations are coupled to a diffusion equation for the 
radiation, the relation between the pressure and den- 
sity is differential. In the Type II instability, it arises 
from the term responsible for generating perpendicular 
radiative fluxes. 

8. The Type II instability describes horizontally propa- 
gating waves that amplify over dynamical time scale. 
If a standing mode is composed out of two oppositely 
propagating modes, it has the appearance of 'bubbles' 
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moving downwards. This arises because the phase ve- 
locity of the wave standing in both the horizontal and 
vertical direction is pointing down. This is opposite 
to the intuitive 'photon bubble' picture, in which the 
structure appears to be moving upwards. 

9. As the flux through an atmosphere increases towards 
the Eddington limit, Thomson scattering dominated 
atmospheres become unstable above a critical fraction 
of the Eddington flux. This instability is independent 
of and comes on top of the possibility of convective in- 
stability. The first mode to become unstable can be 
either Type I or Type II unstable. Type I appears 
at r > 0.5. However, the exact value depends on 
the boundary conditions and thermal gradient. The 
Type II instability appears under more general bound- 
ary conditions but at F w 0.85. 

10. The most unstable horizontal wavelength for the type 
I instability is the longest wavelength in the system at 
Fcrit, but it becomes smaller as F is increased. The 
most unstable mode for the type II instability is about 
0.9 (^Ip^). It doesn't change by much as F is increased 
from Fcrit- 

11. The analysis assumed thus far that the modes can be 
trapped around the layer under consideration. How- 
ever, since the two instabilities operate over a dynami- 
cal time scale, they are probably important even if the 
acoustic oscillations are not trapped, as is the case with 
convection for example. 

12. The unstable modes always tend to reduce the effective 
opacity of the medium (if the luminosity is sufficiently 
close to the Eddington luminosity), thereby increas- 
ing the effective Eddington luminosity. The radiation 
and the opacity must therefore be solved simultane- 
ously and self consistently. Treatment of the inhomo- 
geneities on a large scale (say scales larger than a few 
optical lengths) requires a new formulation of the equa- 
tion of state because of the process of averaging over 
the inhomogeneities. This problem is deferred to a later 
paper. 

All the results presented in this work assume that the 
mass opacity is constant. In general however, the opacity 
can be a function of the gas variables. This can be shown 
to not only change quantitatively the results but also to 
introduce the effects of additional instabilities, some of 
which have not yet been described. Moreover, it should 
be stressed that since the temperature of the matter and 
the radiation at any given point are not necessarily iden- 
tical, non equilibrium effects cannot be neglected. These 
effects should be discussed in the future. 

Finally, if wc arc also to fully understand the behavior 
of the instabilities, we should analyze their nonlinear be- 
havior. It is this non linear behavior that will set many 
of the observational properties of systems shining close to 
their Eddington limit. 

The author wishes to thank the comments and sug- 
gestions made by the editor, anonymous referee, Mitch 
Begelman and Roger Blandford as well as Caltech for the 
DuBridge Prize fellowship which supported him and CITA 
for the current support. 



APPENDIX: A THIRD INSTABILITY? 

The numerical linear analysis performed by Spiegel & 
Tao (1999) (herc^aftc^r S&T) is similar to the analysis car- 
ried out here. However, S&T report on one hand the find- 
ing of a third instability that is not found here. On the 
other hand, they do not find the two instabilities that are 
described in this work. This of course raises the follow- 
ing two questions. Why didn't they find the instabilities 
described here and why didn't the analysis described here 
result with their instability. We wish to address these two 
questions here. 

The answer to the first question is rather simple. 
S&T only analyze one type of boundary conditions. The 
particular boundary conditions they chose tend to pro- 
mote a lowest eigenmode with a net positive correlation 
between the radiation pressure and the gas density. This 
can be seen in their eigenmode profiles by comparing the 
radiation energy density and the gas pressure. Thus, their 
boundary conditions do not allow the 'Type F instability 
to show up. 

The second instability, the 'Type IF, was found to op- 
erate on higher vertical eigenmodes than the lowest one. 
Since S&T used a relaxation method to solve the eigen- 
value problem and not a shooting method as was carried 
out here, they were limited to the lowest eigenmode only. 
Thus, they couldn't find the second instability either. 

The instability found by S&T is similar to various 'over- 
stability' instabilities, such as the K-mechanism, in which 
the growth time scale is significantly longer than the os- 
cillation period. In their case, they found growth rates of 
order lOO's to lOOO's of oscillation periods. Why wasn't 
this instability found here? 

We searched with our code the available parameter 
space and no instability, as the one described, was found. 
One should note that S&T themselves cite the unpublished 
work of Marzek (1977) who studied the same problem but 
could not find this instability either. We could, however, 
artificially get an instability similar to the one described 
in S&T, if we did one of the following: 

1. Calculate the Doppler term wrong, or neglect it alto- 
gether. 

2. Use an inaccurate unperturbed atmosphere. 

3. Have a mismatched opacity between the unperturbed 
atmosphere and the perturbations. 

Since S&T found in their analysis that the Doppler 
terms result with only a very minor correction to the eigen- 
modes, the analysis described in this work was incorrectly 
first executed without the Doppler terms. It was then 
found that the finite speed of light terms - those multi- 
plied by Uc, result with an instability very similar to what 
was described by S&T. When we added the Doppler terms, 
we found that the radiation drag that it induces is always 
comparable but larger than the instability arising from the 
Uc terms. An analytical estimate for the size of the terms 
corroborates this and it shows that the absorption term, 
which introduces damping as well, can be either larger or 
smaller. Namely, there are two possibilities. In the first. 
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S&T calculated the Doppler terms correctly and the origin 
of the instability in not with the Uc terms, in which case 
they must have used a very small Vg/c ratio (found in a 
low temperature atmosphere) a fact which is not stated. 
The second possibility is that the instability does originate 
with the Uc terms in which case they erroneously underes- 
timated the Doppler term. One should note that S&T use 
the Eulerian description of the radiation while we use the 
Lagrangian description, namely, we write the zeroth and 
first moments of the radiation field in the frame of refer- 
ence of the moving material. The former method is, as it 
turns out, significantly more complicated since it involves 
many terms instead of one. 

If the Doppler terms were calculated properly, a second 
way of artificially getting an instability is if the atmosphere 
used as the zeroth order, or unperturbed state, is inaccu- 
rate and doesn't solve the radiation equations accurately. 
S&T describe the exact unperturbed solution, which for 



the optically thin limits is isothermal and the thick limits 
polytropic. They then solve for the eigenmodes assum- 
ing one or the other and then in an atmosphere that has 
both regions. If the latter case is actually not the accu- 
rate solution but the two limits 'glued' together, it would 
explain how an instability can arise. This is because the 
unperturbed state does not solve exactly the unperturbed 
equations. The modes display an instability 'triggered' 
by the atmosphere's will to satisfy the hydrostatic equa- 
tions. We have also found the same eifect for the same 
reasons if the unperturbed atmosphere is calculated with 
only the scattering opacity (which dominates), while the 
linear analysis is calculated with the scattering and very 
small absorption opacity. 

Irrespective of what the resolution to this question is, 
even if the instability does really exist, it is dynamically 
less important than both the Type I and Type II instabil- 
ities. 
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